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1. Survey of theoretical and experimental results in 2D 



Quenched disorder coupled to the energy-density has been the subject of an intensive 
activity in statistical physics, essentially in two-dimensional {2D) systems, during the 
past decades. The qualitative influence of disorder at second order phase transitions 
is well understood since Harris proposed a celebrated relevance criterion [^. At 
first order transitions, randomness obviously softens the transitions, and, under some 
circumstances may even induce second order transitions according to a picture first 
proposed by Imry and Wortis and then stated on more rigorous grounds by Aizenman 
and Wehr P, implying an important result that an infinitesimal disorder induces 
continuous transitions in 2D. The g— state Potts model |^ is the natural candidate for 
the investigations of infiuence of disorder in a perspective linked to critical phenomena, 
since the pure model exhibits two different regimes: a second order phase transition 
when g < 4 and a first order one for g > 4 (in 2D). Many results were obtained in both 
regimes for self-dual quenched randomness in this model in the last ten years, including 
perturbative expansions p, |^, ||, ^ |10|, |lT], |12|, Monte Carlo simulations [|1^, 



|T6|, |T8| or Transfer Matrix calculations [|T9], |0|, |T|, ||, |2|, H, high-temperature 



series expansions |25, 26] or recently short-time dynamic scaling [27|. The numerical 



studies showed that many difficulties, like the lack of self averaging [^, ^ |3y, ^ 
or varying effective exponents due to crossover phenomena may occur. Averaging 
physical quantities over the samples with a poor statistics may thus lead to erroneous 
determinations of the critical exponents. We also note that the previoulsy mentioned 
studies were reported in the case of the random bond system with self-dual probability 
distributions of the coupling strengths in order to preserve the exact knowledge of the 
transition line. 

In real experiments on the other hand, disorder is inherent to the working-out 
process and may result e.g. from the presence of impurities or vacancies in a sample 
produced in Molecular Beam Epitaxy or sputtering experiments. For the description of 
such a disordered system, dilution is thus more realistic than for example a random 
distribution of non-vanishing couplings. Since universality is expected to hold in 
non frustrated random systems, the detailed structure of the Hamiltonian should not 
play any determining role in universal quantities like critical exponents, but dilution 
presumably produces a quite strong disorder compared to e.g. a binary random-bond 
distribution of coupling strengths, and crossover phenomena may alter the determination 
of the universality class. 

Experimentally, the role of disorder in 2D systems has been investigated in 
several systems. Illustrating the inffuence of random defects in the case of the 2D 
Ising model universality class, samples made of thin magnetic amorphous layers of 
(Tbo.27Dyo.73)o.32Feo.68 of 10 A width, separated by non magnetic spacers of 100 A Nb 
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in order to decouple the magnetic layers were produced using sputtering techniques. A 
structural analysis (high resolution transmission electron microscopy and x-ray analysis) 
was performed to characterize the defects inherent to such amorphous structures 
(types of voids or possibly compositional disorder), and in spite of these random 
defects separated on average by a distance of a few nm, the samples were shown to 
exhibit Ising-like singularities with critical exponents (3 = 0.126(20), 7 = 1.75(3) and 
S = 15.1(10) [^. This is coherent with the fact that disorder does not change the 
universality class of the 2D Ising model, apart from logarithmic corrections which are 
not easy to observe experimentally, since their role becomes prominant only in the very 
neighbourhood of the critical point. For this reason, the Ising model is probably not 
the best system to study quantitatively the influence of randomness experimentally. 

More interesting from the point of view of critical phenomena is the case of a 
beautiful experimental confirmation of the Harris criterion reported in a Low Energy 
Electron Diffraction (LEED) investigation of a 2D order disorder transition 
belonging to the 4-state Potts model universality class. Order-disorder transitions of 
adsorbed atomic layers are known to belong to different two-dimensional universality 
classes depending on the type of superstructures in the ordered phase of the adlayer 0, 
p5| . The substrate plays a major role in adatom ordering, as well as the coverage 
(defined as the number of adatoms per surface atom) which determines the possible 
superstructures of the overlayer. For example, sulfur chemisorbed on Ru(OOl) exhibits 
four-state or three-state Potts critical sing ularities for the p(2 x 2) and the ( ^3 x y/3)R30° 
respectively |^ (at coverages 1/4 and 1/2). Other examples can be found in the 
literature, for example oxygen on Ru(OOOl) ordered in a p{2 x 2) |3^ or p{2 x 2)- 
Cs or K on Cu(lll) |^ also belong to the four-state Potts universality class, while 
(\/3 X ^/3)R30° — ICO has three-state Potts singularities The case of the 



(2 X 2)— 2H/Ni(lll) order-disorder transition of hydrogen adsorbed on the (111) surface 
of Ni thus belongs to the 2D four-state Potts model universality class (with expected 
exponents P = 1/12 ~ 0.083, 7 = 7/6 ~ 1.167, and u = 2/3 ^ 0.667 for example). 

Using LEED, it is possible to measure these exponents through the diffracted 
intensity /(q) or structure factor. This is the two-dimensional Fourier transform of 
the pair correlation function of adatom density, where long range fluctuations produce 
an isotropic Lorentzian centered at the superstructure spot position qo with a peak 
intensity given by the susceptibility and a width determined by the inverse correlation 
length, while long range order gives a background proportional to the order parameter 
squared: 



The following exponents were thus measured ^ ^ /9 = 0.11±0.01, 7 = 1.2±0.1 
and z/ = 0.68 ±0.05 in correct agreement with 4-state Potts values (the small deviation. 
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especially for the exponent (3, is possibly attributed to the logarithmic corrections 



to scaling of the pure 4-state Potts model |H2|). The same experiments were then 



reproduced in the presence of intentionally added oxygen impurities, at a temperature 
which is above the ordering temperature of pure oxygen adsorbed on the same substrate. 
The mobility of these oxygen atoms is furthermore considered to be low enough at the 
hydrogen order-disorder transition critical temperature that they essentially represent 
quenched impurities randomly distributed in the hydrogen layer, and the new measured 
exponents become (3 = 0.135 ± 0.010, 7 = 1.68 ± 0.15 and u = 1.03 ± 0.08. It definitely 
rules out the role of these oxygen impurities as extended defects, like steps for example, 
which, according to Ref. |Q, produce a rounding of the transition by a simple finite- 
size-scaling effect. The modification of the universality class in the presence of quenched 
disorder is understood with Harris criterion which predicts such a situation when the 
exponent a of the specific heat is positive for the pure system {a = 2/3 for the 4-state 
Potts model). 

The aim of this paper is to perform numerical simulations of the bond diluted 
Potts model for several values of the number of states per spin (in order to cover the 
two different regimes of the pure system's phase transitions) and provide a reliable 
comparison of the diluted and random-bond problems. This will be achieved through 
a systematic comparison with previous results obtained for self-dual random-bond 
disorder. Here we stress that self-duality is an internal symmetry which can lead to 
some simplifications that the dilute problem does not present and thus it is worth 
comparing both problems. In section ^ the details of the numerical techniques are 
summarized. This section also provides an exposition of the different extrapolation 
techniques in order to leave the discussion of the physical results to the following of the 
paper. Section || presents the phase diagram of the diluted Potts models with g = 3, 4, 
and 8, and section ^jthe results of transfer matrix calculations and the critical behaviour 
of the order parameter correlation function. 



2. Methodology and numerical techniques 

2.1. Definition of the model 

In this paper, we study the 2-dimensional diluted g-state Potts model defined by the 
following Hamiltonian : 

-f3n = Y: K.A.,^, (2) 

where the sum is restricted to nearest neighbours on a square lattice, the degrees 
of freedom {cTj} can take q values and the exchange couplings {Kij} are quenched 
independent random variables chosen according to a binary distribution of non vanishing 
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and vanishing values 

P{K,,) = p6{K,, -K) + {1- p)6{K,,). 
The value p = 1 corresponds to the pure system and Pc 



(3) 



1/2 is the bond 

percolation threshold. Below this threshold there cannot exist any ordered state at 
finite temperature and the phase transition temperature vanishes. 



2.2. Transfer matrix study 

The system is first studied using the transfer matrix method introduced by Blote and 
Nightingale which takes advantage of the Fortuin-Kasteleyn representation 



m 



terms of graphs of the partition function of the Potts model in order to reduce the 
dimension of the Hilbert space. In the Fortuin-Kasteleyn representation, the transfer 
matrix (with no magnetic field) 



Z = Tr + S^^^^Mij), 



(4) 



where Uij = e^'^ — 1, is expanded as a sum over all the possible graphs g (with s sites, 
l{g) loops and c{g) independent clusters) leading to the random cluster model: 



(5) 



bij being the bond variables. Blote and Nightingale suggested to introduce a set of 
connectivity states which contain the information about which sites on a given row 
belong to the same cluster when they are interconnected through a part of the lattice 
previously built. A unique connectivity label rji = rj is attributed to all the sites i of 
such a cluster. In the connectivity space, \Z{m)) is a vector whose components are given 
by the partial partition function Z{m, {f]i}m) of a strip of length m whose connectivity 
on the last row is given by {r^jjm- The connectivity transfer matrix is then defined 
according to \Z{m + 1)) = Tm\Z{m)) and the partition function of a strip of length m 
becomes \Z{m)) = IVk=i '^k\Z{l)), where \Z{1)) is the statistics of uncorrelated spins. 



2.2.1. Free Energy Density The quenched free energy density is given by the 
Lyapunov exponent of the product of an infinite number of transfer matrices 

~ -L-iAo(L). 



fL 

ML) 



lim — In 




ho) 



(6) 
(7) 



where | Vq) is a unit initial vector. In order to determine the error induced by the 
truncation of this product to a finite number of terms, we studied the fluctuations of 
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the free energy density with respect to the number of iterations of the transfer matrix. 
The tests have only been performed for the 4-state Potts model on strips of width L = 8 
with a pseudo critical exchange coupling Kc = 1.16215 at the dilution p = 0.75. This 
value of the exchange coupling is a rough estimate of the critical exchange coupling at 
this dilution as it will be shown later. 




2.06 I ' ' ' ' ' ' ' 1 

250000 500000 750000 1000000 

# iterations 



Figure 1. Convergence of the free energy density (for five different disorder 
realisations) as a function of tlie number of iterations of the transfer matrix. For the 
sake of clarity, error bars are only displayed every 50000 iterations. The parameters 
take the following values L — 8, p = 0.75, Kc = 1.16215. The inset shows statistical 
fluctuations A/ = \J ~ P of free energy density with respect to the number of 
iterations. The dashed line (which is a guide for the eyes) is a power-law with exponent 
-1/2. 

In figure |l], we present the estimates of the free energy density for five independent 
runs and the inset shows its standard deviation. For a self- averaging quantity X, the 
reduced standard deviation squared should be proportional to the inverse volume of the 
system, Rx = (X^— X^)/X^ ~ 1/LN (the number of iterations N of the transfer matrix 
is the length of the strip), and thus at fixed strip width L, AX ~ 1/ \/N . The standard 
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Figure 2. Convergence of the average of the free energy density / obtained by iterating 
10^ times the transfer matrix with respect to the number of disorder reahsations 
entering the average. Each circle corresponds to the free energy density of a given 
sample and the dots to the upper and lower bounds of its error bar. The solid line 
is the running average and the final estimate (/) is written at its right edge. The 
parameters take the following values L = 8, p — 0.75, Kc — 1.16215. 



deviation of the free energy density indeed exhibits an inverse square root decay 1 / ^/N 
with the number of iterations of the transfer matrix (see inset in figure |I|), i.e. it 
reveals that the estimates of the free energy density are weakly correlated from row to 
row and the free energy is self-averaging. It thus turns out to be preferable to average 
the estimates of the free energy density obtained by different runs rather than using very 
long strips for which one would accumulate truncation errors. Figure ^ shows such an 
example of average of the free energy density estimated by independent runs. Despite 
the fact that the standard deviation of the free energy density seems to over-estimate 
the true error by at least one order of magnitude, this definition of the error is a safe 
choice and will be kept in the following. Up to now, all runs have been performed using 
100 independent runs of 10^ iterations of the transfer matrix, leading to an accuracy of 
6-7 digits in the free energy density. 
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2.2.2. Central charge For a pure system, the central charge c is defined as the 
universal coefficient in the lowest-order correction to scaling of the free energy density 

fL. 
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where the regular contribution is 



fo 



(8) 



(9) 



and —u is the exponent associated to the irrelevant vacancy field Eq, 47, 48, 49|. For 



a disordered system, c is defined in the same way from the finite-size behaviour of 
the quenched average free energy density fi, and numerically, since the strip widths 
available are small, we can only expect to measure effective central charges which depend 
on the dilution, Ce//(p), and which would converge towards the true value c in the 
thermodynamic limit. One obviously expects the existence of higher order corrections 
to scaling, but since no analytical expression is known, we cannot include explicite size 
dependence for higher order terms and the possible corrections are taken into account 
by fitting the free energy density including a non universal correction [^, ^ 



Jl — fo — 



nc, 



eff 



+ a2L' 



(10) 



Fits with polynomials in of degrees ranging from 2 to 4 (i.e. up to L ^) were also 
considered, but do not increase the accuracyffl. 



2.2.3. Improvement of the calculation of the central charge In order to obtain more 
accurate estimates of the critical temperature, especially at small values of the number of 
states q where the slow variation of the effective central charge makes this determination 
difficult, we have used a semi- analytical transfer matrix. Instead of computing the 
product of the transfer matrices at a given value of the parameter u = — 1, the 
calculations are performed at u + 6u. All elements of the transfer matrices are expressed 
as polynomials of 6u. This parameter shift 6u is supposed small and the polynomials 

f This choice is arbitrary and has no theoretical grounds, but since we do not know the exact expansion, 
a polynomial fit mimics the effective expansion for the small strip widths available in the numerical 
computations. After several trials, it turns out that the most numerically stable estimates of Ce// are 
those obtained with a polynomial of degree 2 (including only the term) in the range of lattice sizes 
[3:8]. The addition of a term in the same range L e [3 : 8] gives in the vicinity of the maximum 
the same estimate of the central charge within numerical accuracy, i.e. the terms of order higher than 
can be neglected. Polynomials of degree 4 lead to non reliable results, with error bars of order 
Aceff ~ 0.1 due to the small number of degrees of freedom. 
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are truncated to the second order in 6u. The free energy and its error bar are thus 
expressed as polynomials of 6u too and can be numerically computed over a continuous 
range of parameters centered around u. The central charge is then extrapolated using 
the fitting procedure previously presented. An example of such a calculation is given 
by the figure ^ where the solid lines show the central charge (with up and down limits 
of the error bars) as a continuous function of K, compared to a set of independent 
computations (circles) at fixed coupling strengths. 



0.8 



0.79 



0.78 

1.425 1.475 1.525 1.575 

K 



Figure 3. Central charge Ce// of the 3-state Potts model at dilution p = 0.75 versus 
the exchange coupling K as obtained with independent runs of 100 samples (symbols) 
and with a semi- analytic calculation with 500 samples of 10® iterations of the transfer 
matrix (solid curve). The ditted lines correspond to the up and down limits of the 
errors bars of the semi-analytic calculation. 

2.2.4- Phase diagram The phase diagram can be obtained as usually by the evolution 
of the maximum of some diverging quantity, let say the susceptibility for example which 
diverges at criticality in the thermodynamic limit. In the strip geometry considered 
here, we preferred another technique based on the behaviour of the effective central 
charge. In usual critical systems, the Zamolodchikov's c— theorem states that there 
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Figure 4. Dependence of the effective central charge Ce// with the exchange coupling 
K and the bond probability p for the 4-state Potts model. The maximum of Ce// gives 
the location of the transition line. 



exists a c— function decreasing along RG flows and giving the central charge at the fixed 
point In the case of random systems (n in the replica approach), the central 
charge increases and can be expected to reach a maximum value at an optimal disorder 
amplitude. This property, linked to non-unitarity in the presence of disorder, is indeed 
observed in simulations [^, |2^, 0, This property is illustrated in figure || 

where we can follow the maximum of the effective central charge in the plane {p — K) . 



2.2.5. Correlation Functions The spin-spin correlation functions along the strip are 
calculated using an extension of the Hilbert space that allows to keep track of the 
connectivity with a given spin. For a specific disorder realisation, the spin-spin 
correlation function along the strip 

G.(r) = - \ (11) 



where (...) denotes the thermal average, is given by the probability that the spins along 
some row, at columns j and j + r, are in the same state and is expressed in terms of a 
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product of the non-commuting transfer matrices: 



(0 






0) 


(0 




|0) 



(12) 



where | 0) is the ground state eigenvector, T';^, is the transfer matrix in the extended 
Hilbert space. The operator gj identifies the cluster containing aj, while dj+^ gives the 
appropriate weight depending on whether or not crj+r is in the same state as aj. They 
were computed on strips of widths L = 2 to 8 and then averaged over 100 000 disorder 
realisations. 




Figure 5. Histograms of probability distributions of the spin-spin correlation function 
at given distances r — 100, 75, 50 and 25 (from left to right) for the 4-state Potts model 
with an exchange coupling K — 1.16215 and at dilution p = 0.75. Note that the scale 
is logarithmic on the x-axis. 

The numerical calculation of the average over randomness of these spin-spin 
correlation functions is made difficult by the fact that they are not self-averaging. 
Indeed, Figure ^ shows that the probability distribution is close to a log-normal 
distribution (the correlation function is given by a product of matrices whose elements 
are reminiscent of randomness). There are thus some rare events Gaif) ~ (9(1) with a 
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large relative contribution to the average. This is especially true at large distances. An 
accurate estimation of the average requires that the sampling includes such events. As 
can be seen in figure ^ the running average of the spin-spin correlation function presents 
jumps due to these rare events. The error bars are underestimated as the standard 
deviation of the correlation functions (which cannot correspond to the true definition of 
the error on data distributed according to a clearly non Gaussian distribution). 
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Figure 6. Convergence of the spin-spin correlation functions Gair) at a given distance 
r = 100, 75, 50 and 25 (from bottom to top) with respect to the number of disorder 
realisations (up to 6200 samples). 

In the following, results are obtained using an average over 100000 disorder 
realisations. As shown in figure it seems that the rare events have been sufficiently 
well sampled with this number of disorder realisations, since the running average remains 
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Figure 7. Convergence of the average of the spin-spin correlation functions Ga-ir) at 
a given distance r with respect to the number of bins of 6200 disorder realisations. 



flat and no systematic deviation is seen. For some disorder realisations, there can be one 
or several rows with no bond at all. Such cuts disconnect the strip in two subsystems 
and the correlation function is thus vanishing for larger distances for the corresponding 
realisations, inducing a discontinuity of G(j{t). The average being performed over a 
finite number of disorder realisations, these cuts lead to jumps in the average correlation 
functions (figure ^). These jumps are more pronounced for moments G^{t) of increasing 
order n because the contribution of the rare events is enhanced and the average is 
determined by a few configurations that can include such jumps. The averaged moments 
are also less fluctuating (with the distance) for larger strip widths. 
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Figure 8. Jumps in the 4* -moment of the spm-spin correlation functions with respect 
to the distance at q — 4, p — 0.70 and K = 0.81. The different curves correspond to 
lattice sizes ranging from 2 to 8 and the average is performed over 100 000 samples. 
The power L ensures the same asymptotic slope for different strip widths. 



2.3. Methodology used in the computations 

For each value of the number of states, q, a few runs are performed with 100 x 10^ 
iterations in order to have at each strip width an evaluation of the free energy density 
f{K,p). Then, after extrapolation at L ^ oo, the effective central charge Ceff{K,p) 
is evaluated in the temperature-dilution plane. This leads to a first approximate 
determination of the phase diagram Kc{p). For each dilution p, new runs are performed 
with the semi- analytical algorithm (see section 2.2.3|) at Kmax, the coupling strength 



corresponding to the maximum of Cg//. A total of 500 x 10^ iterations of the transfer 
matrix are thus used. When the curve Ce//(p) is very flat, the computation is done for 
two distinct values of p. The spin-spin correlation functions are then computed at the 
maximum of Ceff{p) with 500 000 configurations. 
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3. Phase diagram of the diluted g-state Potts model 

According to the previous section, the phase diagrams of diluted Potts models with 
g = 3, 4 and 8 states per spin are first determined by the location of the maxima of 
the central charge in strip geometries in the {p, K) plane. The effective central charge 
is shown in figure |^ for several dilutions, as a function of K. 
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Figure 9. Estimation of the central charge with respect to the exchange couphng for 
several dilutions (for q = 3 (top), q = 4 (middle) and q = 8 (bottom)). 



The phase diagrams of quenched bond disordered Ising and Potts models were 



studied more than twenty years ago using effective-medium approximation [£4[ [55 
The Hamiltonian is written 

{hi) ihj) 

where Xij = Kij — Km is the deviation from the effective medium homogeneous coupling 
strength. Using the identity Q^^^^'^i'^o = 1 + 5^.^^ (e^*-? — 1), one gets a formal exact 
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Figure 10. Phase diagram of the g-state diluted Potts model obtained numerically 

(symbols) and compared with the single-bond effective medium approximation (solid 
lines). The open circles on each critical line give the loci of the optimal dilutions which 
will be discussed later. 



expression for the thermal average of any quantity Q as 

^ (Qn(,,)[l + ^....(e-'^ -l)])^ . . 

^ ^ (n(,,)[l + 5.,.,(e-^^-l)]>^ ^ ^ 

where {. ■ ■)m stands for the average with Boltzmann factors e~^'^"' with 

-mm=T.KmSa„ay (15) 

(iJ) 

A single bond approximation [xij = everywhere except on one particular bond) then 
leads to the equation of the critical hne 

KM-m^-^^I^fl^S^. (16) 

P-Pc 

This expression should be exact in the vicinity of both the pure system and the 
percolation threshold. Inserting the critical coupling for the pure system Kc — 
ln(l + y/q) = Kc{p — 1) and the percolation threshold Pc — 1/2 of bond percolation 
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on the square lattice indeed leads to an excellent agreement with the numerical data 
(see figure [1^). We also note that the approximation has been improved by a cluster 
extension of the effective interaction [Q. 

4. Critical behaviour 

4.I. Critical behaviour of spin- spin correlation functions 

If the assumption of the existence of a unique stable random fixed point holds, one 
expects that the critical behaviour is asymptotically the same as the system is moved 
along the transition line pc < p < 1- However, in finite systems, one generically has 
to deal with strong crossover effects due to the competition between the disordered 
fixed point and the pure (at p = 1) and percolation (at p = pc) fixed points, or to 
corrections to scaling linked to the appearence of irrelevant scaling variables. It is 
known that these latter effects are generally important in random systems and that the 
corresponding corrections to scaling can be substantially reduced when one measures 
the critical exponents in the regime of the random fixed point, expected to be reached 
at the vicinity of the maximum of the effective central charge in the p-direction (as 
shown below). Let us consider the finite-size behaviour of an observable Q measured 
at a deviation t = K — Kc{p) from the critical point on some system of characteristic 
size L, in the presence of dilution at bond probability p. The variables t and play 
the role of relevant scaling fields (with positive RG eigenvalues yt = l/i' and yi = 1 
respectively), while close to the fixed point, dilution is supposed to be related to some 
irrelevant scaling variable with eigenvalue yp = —uj < 0. At the fixed point there is no 
need that the irrelevant scaling field vanishes, so that one can write p* the corresponding 
dilution and the observable Q obeys the following homogeneity assumption in the scaling 
region 

Q{t,L-\p) = L-^y{L'/%L--{p-p*)) (17) 

which corresponds to the same critical behaviour for any value of p in the range 
Pc < p < 1, described by a unique fixed point. An expansion of the last variable 
(keeping the leading term only) along the critical line (i.e. varying p at Kc{p)) gives 

Q(0, L-\p) = TqL-^^{1 + rg)(p - p*)L-- + ...), (18) 

where the Fq's are non-universal critical amplitudes. It is thus possible to fix p = p* in 
order to minimize the corrections to scaling, and the corresponding value of the dilution 
is empirically found to coincide with the location of the maximum of the central charge 



along the critical line Close to the maximum, the variations of the effective central 
charge itself are small, illustrating that corrections to scaling have there a small influence 
which is consistent with our choice in equation (0). The value of p* is not universal 
and should depend on the system shape, boundary conditions, etc. 
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In figure |Tl|, we show the variation of the effective central charge along the transition 
line. For example in the case g = 4, the random fixed point corresponds roughly to 
the optimum dilution p* ~ 0.700. The estimate of the central charge at this random 
fixed point is c* = 1.0148(40) which is less accurate but perfectly compatible with the 
one obtained in the random-bond case with a binary distribution of coupling strengths 
(c* = 1.0148(4)) |f[ This value has been refined using the semi-analytical transfer matrix 
presented before and a 5-time larger statistic, leading to c* = 1.0143(19) at K* = 1.8072. 
Here and in the following, a star superscript is used to mark the values of the parameters 
which lead to the maximum of Ceff{K,p). 



0.80 



0.79 



q=3 




Figure 11. Estimation of the dilution p*, where the central charge takes its 
maximum value. The data points correspond to computations performed at different 
temperatures, along the critical line, determined by the value of the probability p. The 
solid lines are simple guides for the eyes. 



f We note that a direct comparison between the parameters of dilution (p and Kdp)) with those of 
the random-bond problems is not possible. In the self-dual binary random-bond case for example, two 
couplings K and rK are randomly distributed, leading to a critical line with parameter r, Kdr). 
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Once the optimal values K* andp* are located, we compute the correlation functions 
at the corresponding point in the parameter space. An effective magnetic critical 
exponent Xa{L) for a given strip width L is obtained by fitting the average spin-spin 
correlation functions with the ansatz 

2n 



G,{t) = Aexp (^-— x,(L)rj (19) 

deduced from the logarithmic conformal transformation which maps the infinite plane 
onto an infinitely long cylinder. Here, r is the coordinate along the cylinder (strip 
with periodic boundary conditions in one direction). Conformal mappings of average 
profiles and correlation functions were shown to allow accurate determination of the 
critical exponents in disordered systems, at least in the case of the random-bond Potts 
model [^, and we naturally assume here that such transformations also apply in the 
diluted problem. 

The data for g = 3, 4, and 8 are summarized in table |l|, where our best estimate is 
written in bold face. We can notice that although the critical exponents of the diluted 
and random-bond systems are quite close as expected, the rather small error bars are 
probably underestimated, and may not reflect the true deviation when it is due to 
insufficient disorder average. 

The L-dependence of the magnetic critical exponent is due to corrections to scaling 
(due for example to the above-mentioned crossover effects) which are not explicitly 
taken into account in (p!9|). These effective magnetic critical exponents Xo-(L) are then 
extrapolated in the limit L —>■ +cx3 using a polynomial of degree 2 in 1/L. Strip widths 
between 3 and 8 have been used. The critical exponent associated to the decay of the 
average spin-spin correlation function is for example estimated to be — 0.1419(1) in 
the case g = 4. This estimate is above (but outside error bars) the estimate obtained 
for the 4-state Potts model in the random-bond case (0.1385(3)) in the regime of the 
random fixed point. We also note that a very close value is obtained at the critical 
temperature for the dilution p = 0.725 {xa = 0.1412(1)). 

In the case g = 8, we obviously expect the first-order phase transition of the 
pure 8-state Potts model to be turned into a second order one under the influence of 
randomness [Q]. We indeed observed an exponential decay of the spin-spin correlation 
functions in the longitudinal direction of the strip with a correlation length diverging 
with L, which supports this hypothesis. The phase-diagram is qualitatively the same as 
in the case of the 4-state Potts model. The regime of the random fixed point corresponds 
in this case to the dilution p* ~ 0.65. For the random-bond Potts model, the central 
charge in the regime of the random fixed point was c* = 1.5300(5) which is roughly what 
is obtained in the diluted case at the maximum of the central charge c* = 1.5320(23) at 
dilution p* = 0.650 (figure |Tl|). A calculation using the semi-analytical transfer matrix 
yields the estimate K* = 2.40420 for the critical exchange coupling. Unfortunately, the 
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calculation of the spin-spin correlation functions is made difficult by the large weight of 
rare events leading to decays of correlation functions presenting jumps when the number 
of disorder realisations is too small (as already mentioned at the end of section |^). Using 
500,000 disorder realisations, the decay of the average correlation functions yields the 
critical exponent = 0.1514(2) which is close, although again outside error bars, to 
the estimate in the random-bond case (0.1505(3)). 

Table 1. Central charge and magnetic critical exponent of the g— state diluted Potts 
model compared to the random bond Potts model at the corresponding random fixed 
point. 









Diluted systems 




q 


P 




c 




# of samples 


3 


0.750 
0.775 
0.750 
0.775 


1.5004 
1.4300 
1.50165 

f. 42740 


0.800(4) 
0.799(3) 
0.799(2) 

0.799(2) 


0.13495(6) 

0.13522(6) 


100 X 10*^ 
100 X 10^ 
500 X 10^ 
500 X 10^ 


4 


0.700 
0.725 
0.700 
0.725 


1.8100 
1.7120 
1.8072 
1.7075 


1.0148(40) 
1.0144(40) 
1.0143(19) 
1.0140(18) 


0.1419(1) 

0.1412(1) 


100 X 10^ 
100 X 10^ 
500 X 10*^ 
500 X 10*^ 


8 


0.650 


2.40420 


1.5320(23) 


0.1514(2) 


500 X 10^ 








Random Bond systems 




q 


distribution 


c 






3 


binary** 
ternary 
quaternary** 
continuous** 


0.7998(4) 


0.1347(11) 
0.1344(8) 
0.1343(6) 
0.1344(13) 




4 


binary*^ 




1.0148(4) 


0.1385(3) 




8 


binary*^ 




1.5300(5) 


0.1505(3) 





** From C. Chatelain and B. Berche, Nud. Phys. B 572, 626 (2000). 
From C. Chatelain and B. Berche, Phys. Rev. Em, 3853 (1999). 



In the case g = 3, the procedure is again identical to that of the 4 and 8-state Potts 
model presented before. The peak of the central charge with respect to the exchange 
coupling are narrower than for larger values of q (figure |^). This leads to a better defined 
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phase diagram (figure 0) than for q = 4 and g = 8. Nevertheless, the error bars on the 
central charge are of the order of magnitude of its variation with the dilution p. It is 
thus difficult to define precisely its maximum, approximately located at p* = 0.75. The 
effective critical exponent, as given by the decay of the spin-spin correlation, depends 
much more on the precision on the critical exchange coupling K* than for larger value 
of q. The semi-analytical transfer matrix at the dilution p = 0.75 with an average 
of the free energy over 500 products of 10^ iterations of the transfer matrix gives 
the refined estimate K* ~ 1.50165. The maximum of the central charge is 0.799(2) 
which is compatible, within error bars, with the value 0.7998(4) obtained in the random 
bond case. The estimate of the magnetic critical exponent at this point 0.13495(6) is 
compatible with the value obtained in the random-bond case 0.1347(11) and with the 
estimate 0.13465 obtained by perturbative developments in the neighborhood of g = 2. 



4-2. Multifractal behaviour of the spin-spin correlation function 

In this section, we report a study of the multifractal properties of the spin-spin 
correlation functions of the diluted Potts model. A similar analysis was performed 



in the random-bond case in Ref. p4[. The aim is to provide in the diluted case also 
a test of replica symmetry breaking using the multifractal behaviour of the spin-spin 
correlation functions, then to compare the multifractal spectrum for different values of 

The critical exponent x^-a associated to the algebraic decay of the second moment 



-1/2 

G^{t) of the spin-spin correlation function was calculated using a perturbation 
expansion around the conformal field theory at g = 2 with the assumption that replica 
symmetry holds or is spontaneously broken. At g = 3, the expansion based on the replica 
symmetric scenario was shown to be in a very good agreement with the numerical data 
for the random-bond problem. In the diluted case, we confirm this agreement, as can 
be seen in table 0. 

More generally, the critical behaviour of the moments of the correlation function is 
characterized by a set of exponents which depend on the moment order n in the case 
of multifractality. Numerically, these exponents are obtained by a simple generalization 
of equation (|l^) to: 

-l/n f 271 



Gs(r)' = A„exp(^-— x.n(L)rj (20) 

with the extrapolation to L ^ oo. The exponents associated to the first integer order 
moments are given in table and compared to the corresponding random-bond values 



and to the perturbative results in the replica symmetric scenario [jT3[. The plot is shown 
in figure |12|. 

An interesting quantity is given by the Legendre transform H{a) of the set of 
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Table 2. Magnetic critical exponent of the second moment of the spin-spin correlation 
function of the 3— state diluted Potts model compared to the random bond Potts model 
at the corresponding random fixed point. The perturbative results of Dotsenko et al. 
are recalled. 

Transfer matrix 
Randomness distribution x^i 

Dilution 0.1184(1) 

Random bond binary 0.1177(12)'* 

ternary 0.1182(12)=* 

continuous 0.1173(14)** 

Perturbation 
Replica symmetry 0.11761^ 
Replica symmetry breaking 0.12011^ 

=* From C. Chatelain and B. Berche, Nud. Phys. B 572, 626 (2000). 

^ From Vik. Dotsenko, VI. Dotsenko and M. Picco, Nucl. Phys. B250, 633 (1998). 



Table 3. Magnetic critical exponent of the moments of the spin-spin correlation 
function of the 3— state diluted Potts model compared to the random bond Potts 
model at the corresponding random fixed point. The perturbative results of Lewis are 
recalled (the expansion is only valid for q close to the Ising value q = 2 and at small 
moment order n close to 1). 



Transfer matrix 



n 


theor."* 


binary'' 


continuous'' 


diluted 


1 


0.13465 


0.1347(11) 


0.1344(13) 


0.13522(6) 


2 


0.11761 


0.1177(25) 


0.1173(28) 


0.11841(11) 


3 


0.11006 


0.1051(39) 


0.1070(45) 


0.10552(17) 


4 




0.0938(50) 


0.0996(58) 


0.09511(21) 


5 




0.0822(58) 


0.0906(66) 


0.08630(25) 



** From M.A. Lewis, Europhys. Lett. 43, 189 (1998). 

^ Prom C. Chatelain and B. Berche, Nucl. Phys. B 572, 626 (2000). 
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Figure 12. Exponents of the moments of the correlation function, Xa-^ plotted against 
the moment order for q — 3. A convincing agreement follows from comparison with 
1st and 2nd order expansion calculations and also with numerical results obtained for 
random-bond models. 



exponents X^n = nx^^, 



H [a] = Xfjn, — na, a = — — , n = — - — . [21) 

on da 

This multifractal function is naturally introduced through the probability distribution 
of the quantity of interest. The exponential decay of the moments of the correlation 
function along the strip 

JGA^ = l\G^ir)rPriG)dG ~ exp (^-^X.n(L)r) , (22) 

may be seen as a function Lr{n), which is, by construction, the Laplace transform of 
the probability distribution Prill) at fixed r, with the positive variable y = — lnGcr{T): 



oo 



[G„(r)]" = L,(n) = / e-"^P.(y)dy. (23) 







Inverting the Laplace transform 

-/ 

2i7r Js 



lOO 
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and performing a saddle-point approximation of the Bromwich integral a similar 
exponential expression follows for the probability distribution 

P,(y)~exp(^-^iJ(a)r) (25) 

where the variable a = — 2^ In G'o-(r). We implicitly make use of the assumption 
that the amplitude in equation (|^) only smoothly depends on n||. If the scaling 
dimensions X„n measure the critical decay of the correlation functions (or more generally 
of the moments of the correlation function), the spectral function H{a) measures 
the exponential decay of the correlation function probability distribution, and it is 
absolutely equivalent to work in terms of scaling dimensions or in terms of spectral 
function. The latter quantity is also the natural universal scaling function which 
allows a rescaling of the probability distributions of the correlation functions at different 
distances along the strip or with different strip widths. It is obtained using the identity 
Pr{(y) = (27rr/L)P^(ln Go-(r)) and H{a) = — 2^1nPT-(a), up to an additional term 
which does not depend on a, but does depend explicitly on 2nT/L, and comes from 
the change of variable from y to a and on the possible refinement of the saddle-point 
approximation. 

The maximum of the probability distribution Pr{a) corresponds to the average 
InGo-(r) and defines the typical (or most probable value) of the correlation function. It 
is obtained at the minimum of the spectral function H{a) at position oq. Following the 
definition of a, = l//(2vrr/L) = — lnGo-(T)/(27rr/L) also corresponds to the scaling 
dimension of the typical correlation function. 

27VT 



Gtypir) = exp In G'^(r) = exp — ^"oj • (2^) 

The geometrical interpretation of the function H{a) is the following: the curve H{a) has 
in a a tangent n{x„n — a) of slope —n = OH /da which intercepts the horizontal axis at 
(y. —— Xfj-n and the vertical axis at X At the minimum H{aQ), n = 0, implying 

at the same time that the value of H{ao) also vanishes. For a log-normal distribution, 
the spectral function is simply parabolic, and the deviation from the parabola measures 
the distance from the log-normal probability distribution. 

In the case of the diluted 3— state Potts model, the probability distribution of the 
spin-spin correlation function at distances r = 50, 100, 150, and 200 is presented in 
figure |l^ (in fact Pr{lnG) is shown). As one can see, the distribution is very broad, 
and the broadening is more pronounced at large distances. We mention that the events 
corresponding to vanishing G'o-(t) have been discarded (at r = 50, it corresponds to 3% 
of the events for g = 3 (22% for q = 8) and this proportion increases at r = 200 to 

f a„ being the amplitude of the n*'' order moment, the term — \nan/ {2ttt/ L) vanishes in the hmit of 
large strips t/L ^ 1 and can be forgotten in the inverse Laplace transform. 
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Figure 13. Probability distribution of the spin-spin correlation function for the 3-state 
diluted Potts model at criticality at the optimal dilution p* . a) Probability distribution 
of the In Go-(t) for different distances along the strip, r = 50, 100, 150 and 200 from 
right to left, b) Horizontal rescaling using the variable a = — In Go- (T)/(27rT/i), with 
the same values of r from top to bottom, c) Rough approximation of the universal 
spectral function which should no longer depend on the distance along the strip. It is 
given here at the saddle-point approximation, H{a) — In [2§7^r(cK)] . 



almost 5% for g = 3 (and 52% for g = 8)). A simple change to the natural variable 
a already produces a rescaling in the horizontal direction, and the spectral function 



deduced from the saddle-point approximation, if (a) = — ^^In 



2^P.r(a)J, is shown in 

the same figure. At this approximation, the maxima Pt-(«o) should be given by ^ 
(which take values close 157, 118, 78 and 40 for the values of r chosen here). Since it 
is clearly not true (see figure |T3| b) the minima of if (a) do not vanish and the rescaling 
is not perfect. We thus have to keep quadratic fluctuations about the saddle n* 
leading to the following estimate of the Bromwich integral 



Pr(a) 



exp 



2tit 



H{a] 



(27) 
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0.2 
^ 0.1 



1 1 


1 




Jj - 








Figure 14. Spectral function H{a) for q = 3, 4, and 8 state per spin. The data 
corresponding to distances t = 50, 100, 150 and 200 collapse onto single universal 
spectral functions. The multifractal spectrum becomes wider as q is increased. At 
(7 = 3, the solid lines arc the results deduced from 1st and 2nd order expansions of 
Xct" given in Refs. H and |13|. 



The rescaling is indeed clearly improved with 



H(a) 



2nT 



In 



L 

27TT 



1/2 



provided that the prefactor 
from the maxima of figure 



27r V dn'^ 



-1/2 . 



f28) 



is of order unity. With the numbers taken 
3| b), we deduce that at g = 3 this prefactor takes a value 



close to 1.50, making the correction term to equation (pSf), 2^ In 

-1/2 



2n \ d-n? / 



-1/2 



f For example for g = 3, the prefactors 



In 



take the values 1.48, 1.50, 1.47 



271- \^ dn^ ) 

and 1.44 for r = 50, 100, 150, and 200, respectively, leading to additive corrections at most 

'1/2 ^ 

~ —0.010, which explains the negative shift observed in figure 14 for = 3, 









2tt \ 







while it is even smaller for g = 4 (-0.003) and really negligible for q 
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indeed negligible. This is also the case for g = 4 and 8 and the spectral functions 
presented in figure ^ have a vanishing minimum and produce a very good rescaling of 
the correlation function probability distributions. The multifractal spectrum is broader 
at larger values of the number of states per spin, indicating a more important relative 
weight of the rare events, which implies that numerical results would become less reliable 
for large g's. 

5. Conclusion 

We have studied the magnetic critical properties of the bidimensional diluted Potts 
model for g = 3, 4, and 8. From transfer matrix computations, we obtained the critical 
exponents of the decay of the average spin-spin correlation function, x^, which are close 
to those obtained in the random-bond case. The maximum value of the central charge 
is found to be compatible with the estimates in the corresponding random-bond models. 
Finally, the multifractal spectrum of the diluted 3-state Potts model turns out to be 
compatible with that obtained both perturbatively and numerically for the random- 
bond case. This is a strong evidence of the existence of an unique fixed point for all 
kinds of disorder, even in the absence of duality symmetries. In which concerns the 
multifractal properties, the spectral function H{a) becomes broader when the number 
of states per spin increases, indicating extremely stretched probability distributions of 
the correlation functions. 
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